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ABSTRACT 

We examine the effects of gas-expulsion on initially substructured distributions of 
stars. We perform N-body simulations of the evolution of these distributions in a static 
background potential to mimic the gas. We remove the static potential instantaneously 
to model gas-expulsion. We find that the exact dynamical state of the cluster plays 
a very strong role in affecting a cluster’s survival, especially at early times: they may 
be entirely destroyed or only weakly affected. We show that knowing both detailed 
dynamics and relative star-gas distributions can provide a good estimate of the post¬ 
gas expulsion state of the cluster, but even knowing these is not an absolute way of 
determining the survival or otherwise of the cluster. 
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1 INTRODUCTION 

The vast majority of stars appears to form in environ- 
ments with densities typically much greater than the fiel d 
llLada fc Lada|[20o3 : iBressert et ahllioiol : iKing et ^l2012t) . 
but after a fe w Myr the majority of the stars are dispersed 
into the field (iLada fc Ladall^OSl ). The mechanism by which 
stars are dispersed is unclear, they may form unbound, or 
form in bound clusters which are then unbound by the expul¬ 
sion of residual gas left-over after star formation (see below). 

Star formation does not consume all of the gas in a 
molecular cloud, indeed it is estimate d that at most 30 
per cent of the gas is turned into stars dDobbs et al.l[20l3 : 
IPadoan et al.l l2014l ). Observations show that by 10 Myr, 
and probably well before, young stars are no longer associ¬ 
ated with gas llbada &i LadallioOSi ) . This gas has presumably 
been heated and expelled by feedback (ionisation, mechan¬ 
ical winds and supernovae). It is interesting and important 
to examine how the stars respond to this gas loss and the 
corresponding (very significant?) change in the local poten¬ 
tial. 

Many auth ors have e xamin ed the ef fects of gas loss 
on stars (see iTutukovI lll978t) : iHillsI lll980l ): iMathierJ 
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(1982 

lElmeereenI (Il983ll: iLada. Mareulis Sz DearbornI 

(1984 

): Elmegreen & Glemend (19851): Pinto 

1198 

Z); 

Verschueren & DavidI (19891): Goodwin ( 

1997a 

2 ); 

Gever & Burkerd (200ll): Boilv & Krounal ( 

2003a 

2 ); 

Goodwin & BastianI (200f)l): Bastian & Goodwir 

( 

200 

£); 

Baumgardt & Krounal (20071): Parmentier et al 

( 

200 

= ); 

GoodwirJ (20091)). but the majority of this work has 


concentrated on gas loss from clusters in which the stars 
and gas are both dyn amically relaxed and in glob a l viria l 
equilibrium (but see iLada, Margulis fc DearbornI ll 19841) : 
IVerschueren &: DavidI 1 19891 ): GoodwinI (l2009l )). If one 
assumes that the gas and stars in the cluster are well 
mixed and relaxed then the gas-to-star mass ratio is enough 
to derive the global s tar formation efficien c y (SF E) in 
the clust e r. Ho wever, IVerschueren DavidI (Il989l) and 
I GoodwinI (l2009l ) note that the exact dynamical state of 
clusters at the moment of gas expulsion is also extremely 
important and point-out that the SFE alone is not the most 
important factor in deciding the fate of a cluster. 


Recently it has become very clear that not all 
stars form in relaxed, centrally concentrated clusters, 
and can often form in complex hierarchic al/substructured 


distributions which follow the gas IWhitmore et al. 

1999 

: Ijohnstone et alJ 

2000l: iKirk. Johnstone & Tafalla 

2007 

: iGutermuth et al.l 

I2OO9I: Idi Francesco et al.l 2010l: 
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Schmeia. Kumar fc Ferreira 2008h . It is still unclear which 
of clustered or ‘hierarchical’ is the main mode of star forma¬ 
tion (or how and if they are connected), but it is clear that 
hierarchical is an important mode in many nearby low-mass 
star forming regions. 

In a hierarchical mode of star formation the stars and 
gas are not in anything close to equilibrium. Both will be 
in motion relative to one another as the stars respond only 
to gravity, where as the gas can suffer additional hydrody- 
namical effects and feedback. This means that the absolute 
and relative distributions of stars and gas at formation can 
change significantly before the gas is removed from the sys¬ 
tem. Therefore the effect of gas removal on the stellar distri¬ 
bution depends not just on the relative masses of the stellar 
and gaseous components, but also on how they have dynam¬ 
ically evolved. 

This is the latest in a series of papers in which we have 
examined the resp onse of comp l ex, hierarchic al systems to 
gas expulsion (see ISmith et al.l (l201ll . l2013al B. The com¬ 
plexity of hierarchical distributions means that there is a 
very large parameter space to explore, and a wide variety of 
possible out c omes. In this paper we expand on the work of 
ISmith et al.l ll201ll. l2013al) . 


will be addressed in more detail in a future paper (Blafia et 
ah, in prep.). 


2.1 Initial distributions 

In all cases we model the stellar distribution using N = 1000 
particles with equal masses of 0.5 Mq. 

Using the box fracta l method described by 
iGoodwin fc WhitworthI ll2004l . we create 20 random 
realisations of fractal distributions, each with a fractal 
dimensions of D = 1.6, corresponding to a highly clumpy 
initial distribution within a radius of 1.5 pc. We use the 
same 20 stellar distributions for each background potential. 

We start our simulations with two energies: initial virial 
ratios of Qi = 0.5 (warm), and Qi = 0 (cold). As we will 
show, even our Qi^o.s simulations are not in equilibrium. 
Fractal clusters will then attempt to relax in persuit of equi¬ 
librium and subsequently there are large variations in their 
virial ratio parameter. Thus, we measure Q instantaneously 
at two important epochs: the beginning of the simulation (Qi 
where ‘i’ donates ‘initial’), and the moment when gas expul¬ 
sion begins (Qf where ‘f’ donates ‘final’). The cold systems 
start with the stars initially at rest relative to each-other, 
this is unrealistic, but is the case where we expect the most 
rapid collapse and erasure of substructure. 


2 SIMULATIONS 

Since we wish to continue the work of ISmith et al.l ll201ll . 
1201331) . we use similar, simplified initial conditions for 
our simulation s. We perform N-body simulations using the 
Nbody6 code (lAarsethll2003l) . 

As we describe in more detail below, equal-mass stars 
are distributed in a fractal distribution in a smooth and 
static background potential to mimic the potential of the 
gas in which they are embedded. Given that we use a static 
potential for the gas, we are unable to include active star 
formation in our models. However we do not expect this to 
change the key conclusions of this study. The potential is 
then removed instantaneously to simulate gas-expulsion. 

This is clearly an extreme over-simplification in many 
ways. In reality, the gas is not distributed in a smooth spher¬ 
ical distribution, and both the gas and stars will move in 
response to changes in the global potential. The gas will 
also react to hydrodynamic forces and feedback (which is 
what eventually expel any remaining gas). Gas expulsion is 
unlikely to be instantaneous, rather gas will be lost at dif¬ 
ferent rates in different regions, and dynamics can cause the 
stars and gas to decouple without any feedback. 

We take this simplified approach rather than attempt to 
deal with the gas dynamics with a hydrodynamic method for 
two reasons. Firstly, the practical issue of performing large 
ensembles of simulations - this is much quicker and easier 
with pure A-body simulations. Secondly, the complexity of 
the gas distribution would add large numbers of (largely 
unknown) parameters to our possible parameter space. We 
will return to discuss this issue later. 

We choose equal mass particles in order to avoid com- 
plex two body inte ractions and mass segregation (see e.g. 
lAllison et al.l l2009l for the complications a realistic mass 
spectrum can add to an already complex problem). This 


2.2 The background potential and SFE 

We work with three different static background potentials: 
(i) a Plummer sphere with Rpi = 1.0 pc and Mgas,tot = 
3472 Mq, (ii) a uniform sphere of gas with a maximum ra¬ 
dius of R = 1.8 pc and Mgas.tot = 3455 Mq (equivalent 
to a Plummer sphere with Rpi = cxd ), and finally (hi) a 
highly concentrated Plummer sphere with Rpi = 0.2 and 
Afgas.tot = 2053 Mq . This choice of parameters ensures that 
we obtain a SFE = 0.2 for all three background potentials 
(i.e. we always have exactly 2500 Mq total mass within 
1.5 pc, of which 2000 Mq is gas, and 500 Mq is stars). 

In this work we expel the gas instantaneously at early 
times in the evolution of the cluster i.e. within a few crossing 
times (Iter « 1.4 Myr), and compare to clusters with a later 
gas expulsion (~ 7.5tcr). 

2.3 Gas expulsion time 

We simulate rapid gas expulsion by removing the back¬ 
ground potential instantaneously. T his is the most poten¬ 
tially destructive gas expulsion (see iBaumgardt fc Kroupal 
120071) . 

As we wish to model the effects of early gas expulsion, 
we usually remove the gas potential instantaneously within 
two crossing times {ltd « 1.4 Myr). During this time, 

the initial distributions relax violently and ter may no longer 
be a representative timescale (see section [231l 

We first summarise the results from our previous studies 
before describing and explaining our new results. 

2.4 Previous work 

Numerical models of gas expulsion from initially virialised 
gas-star Plummer spheres have shown that a small fraction 
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of stars can remain bound if the stars make-up more than 
about 30 per cent of the initial system (in this case what 
is often assumed to be a direct measure of the star forma¬ 
tion efficiency), and the majority of the stars will remain 
bound if the fraction is greater than 50 per cent (se e e.g. 


iGoodwin fc BastiarJl2006l : Baumgardt fc Kroupal200'^ ') . The 

speed of gas expulsion is important with fast (instantaneous) 
gas expulsion being significantl y more disruptive than slow 
(adi abatic) gas expulsion (see Baum gardt fc Kroupal l2007l 
and lLada. Margulis fc Dearbort] 19841 '). 

Our initial conditions are a highly simplified, but hope¬ 
fully realistic at a fundamental level, model of a fractal stel¬ 
lar distribution relaxing in a global gas potential before the 
removal of that gas potential. This is very different from the 
initially star-gas equilibrium distributions assumed in most 
previous studies. 

Because the stellar distribution is highly out-of- 
equilibrium and also different from the gas potential, this 
means that the stellar distribution will violently relax. The 
initial fractal substructures will be erased and the stellar 
distribution will become smooth, whilst at the same time 
relaxing to fit the underlying static (g as) potential (see also 
I Allison et al.l[2009l : iParker et al.ll2014i . This means that the 
stellar distribution will become more concentrated relative 
to the static gas potential as potential energy stored in sub¬ 
struc ture is distributed more smoothly (see e.g. lAllison et al] 
| 200 ^ _ 

ISmith et al.l (l201lll identify an important parameter in 
determining the remaining bound fraction as the local stellar 
fraction (LSF). The LSF is a measure of the gas mass within 
the stellar half-mass radius; i.e. a measure of the relative 
importance of the gas to the stars. The LSF is defined as 


LSF = 


M*(r < rh) 


M*(r < rh) -I- Mgas(r < rh) 


( 1 ) 


where rh is the radius that contains half of the total mass in 
stars. M* and Mgas is the mass of stars and gas, respectively, 
measured within rh. 

The LSF is analogous to the star formation efficiency 
(S FE) quoted in many p r evious pape rs (althou g h as noted 
by IVerschueren fc David ill989l ~l and iGoodwinI (l2009l ~l this 
should really be called the effective SFE as its relationship to 
the true SFE is uncertain!. ISmith et al.l (l201ll l show that the 
LSF will depend on the initial distribution of stars, the ini¬ 
tial gas-to-star ma ss, and the initial energy of the stellar dis¬ 
tribu tion (also see IParmentier fc PfalznCTl2013l : IParmentierl 
| 201 ^. _ 

[Smith et al.l (l201ll ') find that there is a reasonably good 
relationship between the final bound fraction and the LSF at 
the point of gas expulsion for systems which ha ve relaxed for 
more t han two initial crossing times. However, ISmith et al.l 
(l2013ah show that if gas expulsion occurs earlier, it is rather 
more complex than this suggests. 

The longer the stars have to relax, the closer to a viri- 
alised, smooth distribution i n equilibrium with th e static gas 
potential they will become. ISmith et al.l (l2013al l show two 
important consequences of this relaxation processes. Firstly, 
the LSF changes with time and so the exact time of gas ex¬ 
pulsion is very important. Secondly, the violent relaxation 
of the initially clumpy stellar distributions is stochastic and 
initial distributions that are initially ‘the same’ (i.e. drawn 
from the same generating functions) can evolve very differ¬ 


ently, and at any particular point in time (i.e. at gas ex¬ 
pulsion) can have quite different dynamics and be at dif¬ 
ferent ‘stages’ in their relaxation. If gas expulsion occurs at 
early times (typically less than one crossing time, or around 
1 Myr) then the LSF ceases to be a good predictor of the 
fina l bound fract i on. 

ISmith et al.l (l2013al ! attempted to quantify these effects 
and found that the virial ratio of the stars at the time of gas 
expulsion is also very important to the final bound fraction 
(as suggested by Goodwin 2009). In this study we concen¬ 
trate on examining the effects of the stellar virial ratio at 
the time of gas expulsion. 


2.5 Motivation 

In thi s paper we extend the work of ISmith et al.l (1201 ll . 
[20133) ■ We have two related questions that we wish to con¬ 
sider. 

Firstly, to what extent is it possible to predict the final 
bound fraction of the system? Secondly, is it ever practically 
possible (either observationally or theoretically) to predict 
the final bound fraction of a particular system? 

In this paper we concentrate mainly on the effects of 
the dynamical state of the stars at the time of gas expulsion 
(as measured by the stellar virial ratio). 

We concentrate on systems which have not had many 
crossing times to relax. For the systems we simulate here 
the instant of gas expulsion are typically 1-3 Myr, or less 
than 2 initial crossing times. This means that the initial 
substructured distributions have not had time to relax and 
are in the process of violent relaxation. It is worth noting 
that this corresponds to the age of the gas-free Orion Nebula 
Gluster (I Jeffries et al.ll201ll L 

We will refer to the virial ratios of the systems, Q — 
T/|II| where T is the kinetic energy, and Q, the potential en¬ 
ergy. Q = 0.5 corresponds to virial equilibrium, but we note 
that our systems (especially initially) are often not in any 
true equilibrium, even if Q = 0.5 (they might be formally 
virialised, but may not have equilibrium spacial or velocity 
distributions). Nevertheless as we shall describe below Q is 
a very useful measure. 


2.5.1 The evolving dynamical state of the cluster 

At the start of the simulation we have a very out-of- 
equilibrium distribution with a Qi = 0 or 0.5. The stars 
will immediately start to violently relax and erase the sub¬ 
structure present i n the system (see also lAllison et al. I I 2 OO 9 I: 
iParker et al.ll20l3! . With our initial conditions there will 
always be an initial collapse of most of the stars. Violent 
relaxation rapidly, but very roughly, attempts to bring the 
system to a rough dynamical equilibrium (both virial equi¬ 
librium of the energies, and a smooth density field). 

Therefore the stellar component of the system rapidly 
changes its density distribution, size, and the way that en¬ 
ergy is distributed. This means that any initial measures 
of size, energy etc. rapidly change, meaning that any useful 
timescale such as crossing time also change. 

We take as a measure of the state of the cluster the 
value of the virial ratio, Q, at any time as well as the rate 
at which the virial ratio is changing, Q. 
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Time [Myr] 

Figure 1. A representative example of the virial ratio variations 
with time of an out-of-equilibrium distribution of stars (a fractal 
in this case) inside a smooth background potential. As we study 
early gas expulsion, the smooth background potential is instan¬ 
taneously removed before two crossing times occur (i.e to the left 
of the vertical green dashed line). 


In Fig. [T] we show the evolution of the virial ratio with 
time for a typical system starting with Qi = 0.5. Even 
though this system starts in ‘virial equilibrium’ it imme¬ 
diately increases its Q, and then oscillates around Q = 0.5 
with decreasing amplitude. 

What happens is that the cluster immediately starts to 
violently relax and attempts to collapse into the gas poten¬ 
tial (thus Q rises as potential energy is converted into ki¬ 
netic energy in the initial collapse). But the initial collapse 
is soon halted and the stellar distribution expands causing Q 
to fall, as the stars oscillate within the potential well of the 
cluster. Whilst this is happening substructure is also being 
disrupted, and within a few oscillations the system smooths 
out and the oscillations represent a ‘pulsation’ of a smooth 
cluster as it attempts to fully virialise. Therefore the oscil¬ 
lations in Q with time provide an internal measure of the 
level to which the system has relaxed. 


2.5.2 Gas expulsion times 

When gas expulsion occurs is (yet another) key parameter 
in setting the final state of the system as quantified by the 
final bound fraction fsee lGoodwirJ[2009l : ISmith et al.H201ll . 
[ 2013 ^. In our simulations this is modelled by the time at 
which we remove the static background gas potential to rep¬ 
resent instantaneous gas expulsion. (Obviously this is a huge 
over-si mplification which w e return to in the discussion.) 

In ISmith et al.l (l2013al l we showed that the value of the 
virial ratio at the start of gas expulsion, Qi, is important - 
is the system in an expandin g or contract i ng par t of its re¬ 
laxation process? However in ISmith et al.l (l2013ah . we chose 
a fixed instant in time for gas expulsion for all fractals. As 
each random realisation of a fractal evolves differently, the 
exact virial ratio at the moment of gas expulsion was very 
varied, and uncontrolled. 

In order to better control the dynamical state of the 
cluster at the point of gas expulsion, we artificially vary the 
instant at which gas expulsion occurs (between 0-2 crossing 
times) so as we can choose the virial ratio of the cluster. 
The upper limit for the time of gas expulsion is marked by 
the green dashed vertical line in Fig.[T] For example, in one 


series of ensembles we ensure that Qi = 0.5 by forcing gas 
expulsion to occur whenever the virial ratio happens to be 
at Q = 0.5. 

Obviously real systems will not always expel gas at a 
pre-determined value of Qf = 0.5, so we also expel the gas 
at other times as Qi varies from subvirial {Qi ~ 0.2) to 
supervirial {Qi ~ 0.7). 


2.6 The full set of initial conditions 

To summarise our set of initial conditions: 

We take ensembles of 10 or 20 statistically identical sys¬ 
tems (all parameters from the same generating functions) 
with N = 1000 equal-mass stars with M = O.hM© dis¬ 
tributed as a D = 1.6 fractal with radius 1.5 pc. The ve¬ 
locities of the stars are scaled to give initial virial ratios for 
the stellar system in the background potential of Qi = 0 or 
Qi = 0.5. 

These stellar distributions sit in a three different static 
background potentials. A Plummer sphere with Rpi = 1 pc, 
a highly concentrated Plummer sphere with Rpi = 0.2 pc, 
and a uniform sphere. All of them ith a total mass of 
2500 M© within 1.5 pc that ensures an effective SFE = 0.2. 

The systems are evolved and their time-evolving virial 
ratios are tracked. The instant of gas expulsion is varied in 
order to have gas expulsion ocur when the final stellar virial 
ratio has a wide range of stellar virial ratios from subvirial 
(Qf ~ 0.2) to supervirial (Qf ~ 0.7). At the moment of gas 
expulsion the local star fraction (LSF) can be calculated. 

They are then evolved until the simulation reaches 15 
Myr (~10.7 initial ter) at which the number of stars still 
bound in a remaining cluster can be found to give the final 
bound fraction, /bound- 

We reiterate that these are not very ‘realistic’ initial 
conditions, and there is much about them that is clearly 
artificial. However, even as artificially simplified as they are, 
they are still extremely messy and complicated. Their use 
is not to model reality directly, but rather to allow us to 
probe the physics behind relaxation and recovery after gas 
expulsion. 


3 RESULTS 


The key parameter that we wish to investigate is the frac¬ 
tion of stars that remain in a bound cluster after gas ex¬ 
pulsion and the post-gas expulsion relaxation of the system. 
This ‘bound fraction’ (/bound) is the size of the naked clus¬ 
ter that remains. To measure the bound fraction we use the 
’snowballing method’. In this technique particles that are 
bound to the cluster are found in an iterative procedure, 
that corrects for the sys temic velocity of th e cluster at each 
iteration, (see sec. 2.6 of ISmith et al. l2013bl for a more com¬ 
plete description). 


3.1 Final bound fractions 

In Figure [2^ we show the final bound fraction, /bound, 
against the local star fraction, LSF, for all the simulations 
we have run in this paper. 

There is a vague trend that a high-LSF results in a high- 
/bound (i.e. the bottom right corner of Figureis empty). 
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Figure 2. a: Crosses show /bound against LSF at the time of gas expulsion for all the simulations carried out in this study (see text 

for details), b: Simulations that are highly sub-virial (Qf = 0.22 -0.24; blue triangles), Qf = exactly 0.5 (green squares), or highly 

super-virial (Qf = 0.68-0.72; red inverted triangles) at the instant of gas expulsion 


But there is a huge amount of scatter in this figure, in partic¬ 
ular around LSF of 0.2 can result in clusters with an /bound 
between zero and almost unity. For any particular value of 
LSF there is a scatter of at least 0.5 in /bound- 

This might suggest that there is no way of estimating 
the final bound fractions of star clusters after gas expulsion. 
We show below that it is possible to understand the system 
and fairly accurately predict the final bound fractions if one 
knows both the LSF and stellar virial ratios at the time of 
gas expulsion. 

Because of how we have (rather artificially) chosen our 
gas expulsion times we can split the simulations shown in 
Figure [2^ into groups depending on their final virial ratios. 
In Figure [5}3 we only plot the simulations with 0.22 < Qf < 
0.24 (blue), Qf ~ 0.5 (green), and 0.68 < Qf < 0.72 (red). 

It is clear from Figure [2 )d that a significant amount of 
the scatter is due to the value of Qf at the time of gas expul¬ 
sion. The 0.22 < Qf < 0.24 simulations all have /bound ~ 1- 
The Qf ~ 0.5 simulations show a rapid increase in /bound 
with LSF for low-LSF, then a very roughly linear increase. 
And the 0.68 < Qf < 0.72 simulations show a roughly linear 
increase in /bound with increasing LSF. 


3.2 A simple physical model 

In Fig. [3] we plot /bound against LSF for bins of different 
Qf increasing from low-Qf in the top left to high-Qf in the 
bottom right. Systems with initial virial ratios of Qi = 0.5 
are marked by filled circles, those with Qi = 0 by open 
circles. 

The black solid lines and blue dashed lines are a simple 
model fit to the data which we describe in this section. Note 
that the colours show the form of the gas potential which 
we will describe in the next subsection. For now we will 
concentrate on building a simple model to fit the /bound 
against LSF trends with different Qf. 

We can construct a very simple analytical model that 
fits the results of our simulations surprisingly well (see 
iBoilv fc Kroupal l2003al for a similar, but rather more de¬ 
tailed derivation). 


As described above and in previous papers, the initial 
fractal stellar distribution will attempt to relax and viri- 
alise within the gas potential. What are important for the 
impact of gas expulsion are two quantities at the time of 
gas expulsion: the virial ratio Qf of the stars relative to the 
gas and the local stellar fraction LSF. The LSF measures 
the relative masses of the gas and the stars within the stel¬ 
lar half-mass radius (see above). Therefore the total mass 
(stars plus gas) Mtot in the region in which the stars are 
present is Mtot ~ M,/LSF. 

Let us denote quantities just before the gas expulsion 
with index 1 and just after the gas expulsion with 2. 

One quantity of interest is the kinetic energy T* of the 
stars, set by their velocity dispersion ct*. If we assume a 
Maxwellian velocity distribution, the kinetic energy is given 
by: 



( 2 ) 


where a is the scale factor of the Maxwellian velocity 
distribution, a is related to the velocity dispersion as a^ = 
(j*7r/(37r — 8). Therefore, 



(3) 


where n = 7r/(37r-l-8). After gas expulsion the stars have 
not had time to change their kinetic energy (since the gas is 
expelled instantaneously) and so we can assume r *,2 = 7»,i- 
The potential energy of the stars before gas expulsion 
can be approximated by 


fl*,i 


-M. 


GMtot 

Fh 


(4) 


were G is Newton’s gravitational constant, while the poten¬ 
tial energy after the gas is lost is only due to the potential 
made by the stars alone 


11*2 ~ -M*^^ = LSFfI*i. 

Fh 


(5) 


Now we calculate the escape velocity of the system after 
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Figure 3. The /bound"LSF trend for different virial ratios. Colors represent the shape of the background gas, a Plummer sphere (blue) 
and a uniform sphere (red), filled circles are simulations with Q\ = 0.5 and open circles are distributions with Q\ = 0.0. The black solid 
lines and blue dashed lines are the fits from the model described in sec. 13.21 


the gas is gone as 



If we now replace by LSF times 
have 


( 6 ) 

-T*,\IQi we 


t’esc 



(7) 


A reasonably first guess of /bound would be the fraction 
of stars with velocities below the escape velocity. If we as¬ 
sume a Maxwellian velocity distribution, then /bound is given 
by its cumulative probability distribution with the form: 


F«.V) = erf(^x)-^Xexp(-^) (8) 

where X = Vescja. Since = kct* then X = Uesc/%/«'’'* 
and finally: 


/bound 


erf 


f /3LSF^ /6 LSF 
[y 2 Qs ) y TV Qt 


exp 


3LSF=\ 

2g? J 


(9) 


In Fig. [3 we show /bound against LSF for various values 
of Q{. The solid black line is the fit from above which has 
no free parameters. This simple model describes the data 
points of our simulations very well, especially if we look at 
high LSF and low Qf values, i.e. when we do not lose many 
unbound stars (upper panels). 


When we have high Qf values as in the lower panels 
of Fig. [3] the simple model (solid black line) tends to over¬ 
estimate the final bound fraction. We can apply a simple 
correction. Following the first estimate of /bound a fraction 
of stars is lost very rapidly after gas expulsion, and so the 
escape velocity falls by a further factor V /bound in Eq. [3 We 
then have to solve Eqs. [3 and |8] iteratively which gives the 
blue dashed-lines in Fig. [3l In most cases the true values of 
/bound are enclosed between the solid black and blue dashed- 
lines suggesting that reality is somewhere inbetween. 

We have constructed a simple analytic approximation 
with no free parameters that estimates the final bound frac¬ 
tion from the values of the stellar virial ratio and LSF at the 
moment of gas expulsion. Given the simplifying assumptions 
we have made it is very gratifying that this seems to explain 
the results so well. 


3.3 The effect of the gas potential 

In Fig.|3points are coloured according to the form of the gas 
potential: blue is a Plummer potential, and red a uniform 
sphere. There appears to be a very strong dependency on 
the form of the gas potential. In Fig. 3 systems with con¬ 
centrated gas potentials {Rpi = 1 pc) shown by the blue 
markers are concentrated to the left of each panel with low 
LSF and low /bound. Systems with extended gas potentials 
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{Rpi = oo) shown by the red markers are towards the right 
with higher LSF and /bound- 

Taken at face value this suggests that the form of the 
gas potential is crucial in determining the fate of a system. 
However, this is not the case. Rather it is due to a link be¬ 
tween the form of the gas potential and the possible values 
of the LSF. The LSF measures the relative masses of gas 
and stars within the half-mass radius of the stars. Gas out¬ 
side this radius is not taken into account. Even though the 
total mass in gas in the whole star forming region stays con¬ 
stant, the LSF fluctuates as the half-mass radius of the stars 
fluctuat es (this is th e moti vation for the introduction of the 
LSF by ISmith et akl (1201 il l. 

In a bound, fractal distribution the stars can do nothing 
except collapse to a denser (and smoother) configuration. 
Much of the initial potential energy in a fractal distribu¬ 
tion is localised in substructure and is redistributed during 
violent relaxation. The potential energy, H, of a system is 


Q 


-A 


GM^ 

R 


( 10 ) 


where M is the mass of the system and R some charac¬ 
teristic radius (and G the gravitational constant). A is a 
measure of the mass distribution of the system. For a Plum¬ 
mer sphere, if R is the Plummer radius then A ~ 0.3. But 
for a, D = 1.6 fractal, when R is the initial size of the sys¬ 
tem, A ~ 1.5. Therefore, the violent relaxation of a fractal 
causes a significan t decrease in the size of the system (see 
lAllison et al.|[200^ for details). 

Exactly how such a system will contract depends on 
the exact details of the initial fractal distribution, the initial 
virial ratio (Qi = 0 systems will contract more than Qi = 0.5 
systems) and how relaxed the system has become. However, 
we Hnd it does not depend on the shape of the background 
gas potential, as shown in Fig. [I] In the upper panel, symbols 
with error bars are the average LSF of simulations at the 
moment of gas expulsion. We include data points for (from 
left to right) Rpi = oo (uniform gas), Rpi = 1 pc, and also 
Rpi — 0.2 pc (a very concentrated gas distribution). On the 
x-axis, we plot 1/Rpi in order to place all simulations on 
the same plot. There are two curves for the two different 
initial virial ratios {Qi = 0.0 and 0.5). There is a clear trend 
for the LSF to be lower as the gas becomes increasingly 
concentrated. To understand this, we must bear in mind that 
the LSF is a function of the total gas mass within the half 
mass radius of the stars. Therefore a change in LSF could 
arise from either a change in the amount of gas surrounding 
the stars, and/or a change in the half-mass radius of the stars 
as the gas scale length is varied. We find that the half-mass 
radius is only a very weak function of the gas scalelength as 
shown in the lower panel of Fig. [d] Here symbols with error 
bars are the average half-mass radius Rh of the stars. This 
weak dependency demonstrates that the strong dependency 
of the LSF on gas scalelength arises mainly for the following 
reason - by making the gas more concentrated, more gas is 
being placed about the stars, and the LSF is lowered. 

To conhrm that the small variation in Rh with gas scale- 
length does not play a strong role, we calculate the average 
Rh for each set (see horizontal dashed lines in the bottom 
panel of Fig. H- Now we fix Rh to have the average value 
(i.e a constant value for all gas scalelengths) and recalculate 
the LSF values at their new half-mass radius. The results 



1 / [pc '] 


Figure 4. The variation of the LSF of the clusters due to the 
change in the concentration of the background gas. Top panel: 
The average of the LSF of the simulations against the inverse of 
their scale lengths. A black solid line connects simulations with 
the same initial virial ratio as labelled. Bottom panel: The half 
mass radius is only weakly dependent on the gas scale length. 
The average is shown by the horizontal dashed line. The brown 
dashed line in the upper panel is the recalculated LSF using a 
fixed half-mass radius with the average value. 


are indicated by the brown dashed lines in the upper panel. 
The trend of LSF with gas scalelength is very similar, even 
when the stellar half-mass radius is fixed to be constant. 
This confirms that the strong dependency of the LSF on gas 
scalelength arises almost entirely for the following reason. 
Increasing the gas concentration places more gas about the 
stars, and does not change the stellar distribution signifi¬ 
cantly. 


4 DISCUSSION AND CONCLUSIONS 

Initially clumpy and irregular distributions of stars cannot 
be in dynamical equilibrium. As a result, they undergo vi¬ 
olent relaxation with initially significant changes in their 
virial ratio as they expand and collapse, attempting to ap¬ 
proach equilibrium. This occurs even when the clusters are 
initially ‘virialised’ (ie. Qi = 0.5). These deviations are 
largest for very young star clusters, and decrease as the clus¬ 
ter settles down, as substructure is erased. As a result, the 
effects of gas expulsion at early times, before the system has 
relaxed, depend strongly on the instantaneous value of the 
virial ratio as well as the Local Star Fraction (LSF, relative 
distribution of stars in the gas potential). 

At later stages (>2 crossing-times), it is known that 
the LSF becomes the key predictor of cluster survival from 
gas expulsion, with secon d-order modification s due to the 
cluster’s dynamical state (ISmith et al.ll2013al l. However at 
these early stages when oscillations in the virial ratio are so 
large, we have shown that the dynamical state of the cluster 
may actually be equally influential (if not more influential) 
than the LSF. 

A primary goal of studying the response of young, em¬ 
bedded star clusters to gas expulsion is to predict how well 
a cluster survives gas expulsion, based on its pre-gas expul- 
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sion properties. This study reveals that both the LSF, and 
the dynamical state can be important parameters dictating 
cluster survival to gas expulsion. Fortunately in our numeri¬ 
cal studies, we can ascertain the exact value of the LSF and 
virial state. However, observationally, it may be incredibly 
challenging to measure either of these properties accurately. 
It is not inconceivable that the LSF might be calculated ap¬ 
proximately by deprojection, although it would need to be 
a cluster caught very close to the instant of gas expulsion, 
or the LSF may later change. However, measuring the virial 
ratio of a real cluster is a huge challenge. 

To worsen matters, our study reveals that in certain 
circumstances, even with a knowledge of both the LSF and 
virial ratio, the cluster survival maybe poorly constrained. 
For example, take a cluster with a low virial ratio (e.g 
(5f=0.34 at gas expulsion; upper-left panel of Fig. |3}. If the 
cluster has an LSF~0.2 (a reasonable physical value), the 
/bound-LSF trend rises very steeply. Such a cluster is equally 
likely to be near destroyed (have ~ 90% of its stars un¬ 
bound), as only weakly affected (losing ~ 30% of its bound 
stars). Thus it is possible that, even if the virial ratio were 
measured, the result could place the cluster in a region of pa¬ 
rameter space where the cluster survival could be anything 
from weak mass loss to near total destruction. 

Comparing the panels of Fig. [21 we can see that clusters 
with LSF~0.2-0.4 are the most sensitive to their dynamical 
state. In comparison clusters with high LSF vary their re¬ 
sulting bound fractions very little, even for large changes in 
dynamical state. If LSF~0.2 is a typical value, then these 
results suggest that clusters which are observed post-gas ex¬ 
pulsion, must have been sub-virial to avoid losing a large 
fraction of their stellar mass during the process. 

Clearly our models are extremely simple conception- 
ally. They lack a large number of physical processes that 
are also highly important in young star clusters. For exam¬ 
ple, our cluster stars have no initial mass function, we start 
our simulations with no binaries, we do not consider stellar 
evolution, and our treatment of the gas is highly simpli¬ 
fied. Nevertheless, the use of such simple idealised models 
has enabled us to clearly determine the significant role of 
clumpy substructure and the dynamical state of the clusters 
on cluster survival following gas expulsion, through the use 
of controlled numerical experiments. This approach has re¬ 
vealed just how sensitive star clusters are to their dynamical 
state when gas expulsion occurs. We therefore suggest that 
real star clusters will be very sensitive, perhaps as sensitive 
as our model star clusters, to their dynamical state when 
the gas is expelled at early times. 

Our key results may be summarised in the following: 

(i) For early gas expulsion (before 2 crossing-times) we 
find the dynamical state of our model star clusters, measured 
at the time of gas expulsion, plays a key role in influencing 
cluster survival following gas expulsion. Star clusters may 
be highly super- or sub-virial in these early phases. 

(ii) We show how the /bound-LSF trend can be well ap¬ 
proximated with a very simple analytical model. The model 
matches the simulations best when the dynamical state is 
not extreme (i.e highly super- or sub-virial). 

(iii) Clusters which have LSFs in the range 0.2-0.4 (phys¬ 
ically reasonable values) are most sensitive to the virial ratio 
at the instant of gas expulsion. 


(iv) Clusters with low virial ratio have a very steep rise in 
the /bound-LSF trend. For such a cluster with an LSF~ 0.2, 
it is therefore not possible to predict if the cluster will be 
heavily destroyed or only mildly affected - even knowing 
both the LSF and the virial ratio. 

This study highlights the difficulties faced in trying to 
determine the survival rate(s) of real star clusters due to gas 
expulsion. At early times, the dynamical state of a cluster 
may be far from dynamical equilibrium, and this can signif¬ 
icantly affect the clusters survival to gas expulsion. Thus a 
best estimate of a cluster’s survival is found measuring both 
the LSF and virial ratio. Accurately measuring these two 
parameters for a real cluster represents a huge observational 
challenge, in particular the dynamical state. Furthermore, 
some clusters may be situated in regions of parameter space 
where their survival to gas expulsion remains highly uncer¬ 
tain, even knowing both the LSF and virial ratio. 
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